Compiled under R version 3.6.3 (2020-02-29)
WARNING: edit the working directory to your preferred folder.
This document details all analyses performed in R for the study:
Legendre, L. J., D. Rubilar-Rogers, G. M. Musser, S. N. Davis, R. A. Otero, A. O. Vargas, and J. A. Clarke. 2020. A giant soft-shelled egg from the Late Cretaceous of Antarctica. Nature 583, 411–414.
For more information regarding the study, datasets, and analyses, please refer to the Supplementary Information of this paper. If you have any additional questions, feel free to email me at lucasjlegendre@gmail.com.
library(AICcmodavg)
library(ape)
library(caper)
library(dplyr)
library(evobiR)
library(ggplot2)
library(MPSEM)
library(nlme)
library(nortest)
library(phylopath)
library(phytools)
library(RColorBrewer)
This part requires the use of dataset 1 (“Dataset1-lepidosaurs.txt”; see also Supplementary Table 1 in Legendre et al., 2020), and lepidosaur phylogenetic tree (“Lepidosaurtree.trees.nex”).
tree<-read.nexus("Lepidosaurtree.trees.nex")
data2<-read.table("Dataset1-lepidosaurs.txt", header=T); data2V<-subset(data2, !is.na(V)&!is.na(SVL))
pruned.tree<-drop.tip(tree, setdiff(tree$tip.label, data2V$Species))
data<-read.table("Dataset1-lepidosaurs.txt", header=T, row.names = "Species"); dataV<-subset(data, !is.na(V)&!is.na(SVL))
alpha <- seq(0, 1, 0.1)
fit <- list()
form<-log(V)~log(SVL)
for (i in seq_along(alpha)) {
cor <- corMartins(alpha[i], phy = pruned.tree, fixed = T)
fit[[i]] <- gls(form, correlation = cor, data = data, na.action=na.exclude, method = "ML")
}
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
sapply(fit, logLik)
## [1] Inf -325.1206 -325.9100 -326.0731 -326.1266 -326.1503 -326.1665
## [8] -326.1806 -326.1939 -326.2064 -326.2182
g <- seq(0.1, 1, 0.1)
fit <- list()
form <- log(V)~log(SVL)
for (i in seq_along(g)) {
cor <- corBlomberg(g[i], phy = pruned.tree, fixed = T)
fit[[i]] <- gls(form, correlation = cor, data = data, na.action=na.exclude, method = "ML")
}
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
sapply(fit, logLik)
## [1] -322.9247 -331.3311 -337.6858 -341.8988 -344.8120 -346.9305 -348.5380
## [8] -349.8001 -350.8185 -351.6589
BM<-gls(log(V)~log(SVL), data=dataV, correlation=corBrownian(phy=pruned.tree), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OU<-gls(log(V)~log(SVL), data=dataV, correlation=corMartins(0.1, phy=pruned.tree, fixed=T), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
Lambda<-gls(log(V)~log(SVL), data=dataV, correlation=corPagel(1, phy=pruned.tree), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
EB<-gls(log(V)~log(SVL), data=dataV, correlation=corBlomberg(0.1, phy=pruned.tree, fixed=T), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OLS<-gls(log(V)~log(SVL), data=dataV, method="ML")
Cand.models = list()
Cand.models[[1]] = BM
Cand.models[[2]] = OU
Cand.models[[3]] = Lambda
Cand.models[[4]] = EB
Cand.models[[5]] = OLS
Modnames = paste(c("BM", "OU", "Lambda", "EB", "OLS"), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = T)
summary(Lambda)
## Generalized least squares fit by maximum likelihood
## Model: log(V) ~ log(SVL)
## Data: dataV
## AIC BIC logLik
## 641.6771 655.6163 -316.8385
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.5664802
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1.978195 0.5517311 3.585434 4e-04
## log(SVL) 1.114992 0.0559233 19.937869 0e+00
##
## Correlation:
## (Intr)
## log(SVL) -0.51
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.0709107 -0.7544921 -0.1511364 0.4256780 2.0898411
##
## Residual standard error: 1.157687
## Degrees of freedom: 241 total; 239 residual
plot(Lambda, resid(., type="n")~fitted(.), main="Normalized Residuals v Fitted Values",abline=c(0,0))
res<-resid(Lambda, type="n")
par(mar=c(5.1,4.1,4.1,2.1))
qqnorm(res)
qqline(res)
lillie.test(residuals(Lambda))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuals(Lambda)
## D = 0.053479, p-value = 0.09167
lillie.test(chol(solve(vcv(pruned.tree)))%*%residuals(Lambda))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: chol(solve(vcv(pruned.tree))) %*% residuals(Lambda)
## D = 0.078349, p-value = 0.001089
which(res>4); which(res<(-2.5))
## named integer(0)
## Chirindia_ewerbecki Lycodon_striatus Micrurus_mipartitus
## 89 196 222
## Ogmodon_vitianus Leptotyphlops_scutifrons Myriopholis_blanfordi
## 224 227 228
## Indotyphlops_braminus
## 235
newdata<-dataV[-c(89,168,222),]
newdata2<-data2V[-c(89,168,222),]
newtree<-drop.tip(pruned.tree, setdiff(pruned.tree$tip.label, newdata2$Species))
Newlambda<-gls(log(V)~log(SVL), data=newdata, correlation=corPagel(1, phy=newtree), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
summary(Newlambda)
## Generalized least squares fit by maximum likelihood
## Model: log(V) ~ log(SVL)
## Data: newdata
## AIC BIC logLik
## 615.7252 629.6143 -303.8626
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.3749586
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 1.773470 0.4386847 4.042699 1e-04
## log(SVL) 1.137279 0.0542297 20.971523 0e+00
##
## Correlation:
## (Intr)
## log(SVL) -0.618
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -3.4479913 -0.7458697 -0.0785325 0.6087246 2.5045329
##
## Residual standard error: 0.9939715
## Degrees of freedom: 238 total; 236 residual
plot(Newlambda, resid(., type="n")~fitted(.), main="Normalized Residuals v Fitted Values",abline=c(0,0))
res<-resid(Newlambda, type="n")
par(mar=c(5.1,4.1,4.1,2.1))
qqnorm(res)
qqline(res)
lillie.test(residuals(Newlambda))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuals(Newlambda)
## D = 0.057905, p-value = 0.05125
lillie.test(chol(solve(vcv(newtree)))%*%residuals(Newlambda))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: chol(solve(vcv(newtree))) %*% residuals(Newlambda)
## D = 0.091664, p-value = 4.941e-05
caper#newlambda:
datacomp<-comparative.data(phy=newtree, data=newdata2, names.col="Species")
pgls<-pgls(log(V)~log(SVL), data=datacomp, lambda="ML")
summary(pgls)
##
## Call:
## pgls(formula = log(V) ~ log(SVL), data = datacomp, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.107910 -0.039203 -0.005097 0.034535 0.147102
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.718
## lower bound : 0.000, p = 1.8499e-05
## upper bound : 1.000, p = 4.5223e-05
## 95.0% CI : (0.413, 0.908)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.41450 0.91690 -0.4521 0.6528
## log(SVL) 1.63215 0.15212 10.7293 6.661e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06154 on 64 degrees of freedom
## Multiple R-squared: 0.6427, Adjusted R-squared: 0.6371
## F-statistic: 115.1 on 1 and 64 DF, p-value: 6.121e-16
#lambda:
datacompl<-comparative.data(phy=pruned.tree, data=data2V, names.col="Species")
pgls<-pgls(log(V)~log(SVL), data=datacompl, lambda="ML")
summary(pgls)
##
## Call:
## pgls(formula = log(V) ~ log(SVL), data = datacompl, lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.140589 -0.048249 -0.004924 0.031043 0.155446
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.715
## lower bound : 0.000, p = 7.7636e-05
## upper bound : 1.000, p = 0.00015352
## 95.0% CI : (0.381, 0.918)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.19647 0.96534 0.2035 0.8393
## log(SVL) 1.53068 0.16093 9.5115 5.418e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.06683 on 66 degrees of freedom
## Multiple R-squared: 0.5782, Adjusted R-squared: 0.5718
## F-statistic: 90.47 on 1 and 66 DF, p-value: 5.417e-14
Not much difference between the two models, assumptions of normality are not met in either of them – we use Lambda instead of Newlambda for the plot.
ggplot(dataV, aes(log(SVL), log(V), color=Suborder)) +
geom_point(size=4) +
xlab("ln snout-vent length (mm)") +
ylab("ln egg volume (mm3)") +
geom_abline(intercept=Lambda$coefficients[1], slope=Lambda$coefficients[2], colour="lightblue", size=1.3) +
theme(panel.background = element_rect(fill="black")) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_colour_brewer("Suborder", palette="Accent")
exp((-0.7990941/1.6708400)+(1/1.6708400)*log(5471405.789))
## [1] 6684.303
data[1,7]<-6684.303
dataV2<-subset(data, !is.na(V))
ggplot(dataV2, aes(log(SVL), log(V), color=Suborder)) +
geom_point(size=4) +
xlab("ln snout-vent length (mm)") +
ylab("ln egg volume (mm3)") +
geom_abline(intercept=Lambda$coefficients[1], slope=Lambda$coefficients[2], colour="lightblue") +
theme(panel.background = element_rect(fill="black")) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_colour_brewer("Suborder", palette="Accent")
Eggshell thickness ~ body mass is not detailed here, since the code is virtually identical – only variable names and regression coefficients for the plot and estimation are different.
LambdaBM<-gls(log(V)~log(BM), data=dataV, correlation=corPagel(1, phy=pruned.tree), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
summary(LambdaBM)
## Generalized least squares fit by maximum likelihood
## Model: log(V) ~ log(BM)
## Data: dataV
## AIC BIC logLik
## 524.34 538.2792 -258.17
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.3112388
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 5.260359 0.26352407 19.96159 0
## log(BM) 0.645885 0.02268114 28.47677 0
##
## Correlation:
## (Intr)
## log(BM) -0.297
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.91042516 -0.83939318 -0.05097508 0.63224435 3.05176705
##
## Residual standard error: 0.7847368
## Degrees of freedom: 241 total; 239 residual
exp(-(5.050557/0.619536)+(1/0.619536)*log(5471405.789))
## [1] 21657214
data[1,8]<-21657214
dataV2<-subset(data, !is.na(V))
ggplot(dataV2, aes(log(BM), log(V), color=Suborder)) +
geom_point(size=4) +
xlab("ln body mass (g)") +
ylab("ln egg volume (mm3)") +
geom_abline(intercept=LambdaBM$coefficients[1], slope=LambdaBM$coefficients[2], colour="lightblue") +
theme(panel.background = element_rect(fill="black")) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_colour_brewer("Suborder", palette="Set1")
dataVT2<-subset(data2, !is.na(Thickness)&!is.na(V)); dataVT2<-dataVT2[-1,]
pruned.treeVT<-drop.tip(tree, setdiff(tree$tip.label, dataVT2$Species))
dataVT<-subset(data, !is.na(Thickness)&!is.na(V)); dataVT<-dataVT[-1,]
BMVT<-gls(log(Thickness)~log(V), data=dataVT, correlation=corBrownian(phy=pruned.treeVT), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OUVT<-gls(log(Thickness)~log(V), data=dataVT, correlation=corMartins(1, phy=pruned.treeVT), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
LambdaVT<-gls(log(Thickness)~log(V), data=dataVT, correlation=corPagel(1, phy=pruned.treeVT), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Error in corFactor.corStruct(object): NA/NaN/Inf in foreign function call (arg 1)
EBVT<-gls(log(Thickness)~log(V), data=dataVT, correlation=corBlomberg(0.1, phy=pruned.treeVT, fixed=T), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OLSVT<-gls(log(Thickness)~log(V), data=dataVT, method="ML")
Cand.models = list()
Cand.models[[1]] = BMVT
Cand.models[[2]] = OUVT
Cand.models[[3]] = LambdaVT
## Error in eval(expr, envir, enclos): object 'LambdaVT' not found
Cand.models[[4]] = EBVT
Cand.models[[5]] = OLSVT
Modnames = paste(c("BMVT", "OUVT", "LambdaVT", "EBVT", "OLSVT"), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = T)
## Error in formatCands(cand.set):
## Functions do not support mixture of model classes
summary(EBVT)
## Generalized least squares fit by maximum likelihood
## Model: log(Thickness) ~ log(V)
## Data: dataVT
## AIC BIC logLik
## 147.5644 154.2229 -70.78221
##
## Correlation Structure: corBlomberg
## Formula: ~1
## Parameter estimate(s):
## g
## 0.1
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 2.3502382 0.4486431 5.238548 0
## log(V) 0.2840476 0.0555982 5.108931 0
##
## Correlation:
## (Intr)
## log(V) -0.966
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.6576055 -0.5886486 -0.1688905 0.5376755 2.7915997
##
## Residual standard error: 0.7289907
## Degrees of freedom: 68 total; 66 residual
lillie.test(residuals(EBVT))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuals(EBVT)
## D = 0.071519, p-value = 0.5269
lillie.test(chol(solve(vcv(pruned.treeVT)))%*%residuals(EBVT))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: chol(solve(vcv(pruned.treeVT))) %*% residuals(EBVT)
## D = 0.10736, p-value = 0.05007
plot(EBVT, resid(., type="n")~fitted(.), main="Normalized Residuals v Fitted Values",abline=c(0,0))
res<-resid(EBVT, type="n")
par(mar=c(5.1,4.1,4.1,2.1))
qqnorm(res)
qqline(res)
which(res<(-2)); which(res>2.5)
## Zootoca_vivipara_viviparous Crotalus_viridis
## 45 67
## Chamaeleo_senegalensis
## 19
newdataVT<-dataVT[-c(19,45),]
newdataVT2<-dataVT2[-c(19,45),]
newtreeVT<-drop.tip(pruned.treeVT, c("Zootoca_vivipara_viviparous","Chamaeleo_senegalensis"))
NewEBVT<-gls(log(Thickness)~log(V), data=newdataVT, correlation=corBlomberg(0.1, phy=newtreeVT, fixed=T), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
summary(NewEBVT)
## Generalized least squares fit by maximum likelihood
## Model: log(Thickness) ~ log(V)
## Data: newdataVT
## AIC BIC logLik
## 126.345 132.9139 -60.17249
##
## Correlation Structure: corBlomberg
## Formula: ~1
## Parameter estimate(s):
## g
## 0.1
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 2.343769 0.3932038 5.960697 0
## log(V) 0.284315 0.0482883 5.887863 0
##
## Correlation:
## (Intr)
## log(V) -0.965
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.4180308 -0.6450608 -0.1851669 0.6097965 2.0464452
##
## Residual standard error: 0.6409687
## Degrees of freedom: 66 total; 64 residual
plot(NewEBVT, resid(., type="n")~fitted(.), main="Normalized Residuals v Fitted Values",abline=c(0,0))
res<-resid(NewEBVT, type="n")
par(mar=c(5.1,4.1,4.1,2.1))
qqnorm(res)
qqline(res)
lillie.test(residuals(NewEBVT))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: residuals(NewEBVT)
## D = 0.070589, p-value = 0.5719
lillie.test(chol(solve(vcv(newtreeVT)))%*%residuals(NewEBVT))
##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: chol(solve(vcv(newtreeVT))) %*% residuals(NewEBVT)
## D = 0.088221, p-value = 0.2293
NewEBVT$coefficients
## (Intercept) log(V)
## 2.343769 0.284315
ggplot(newdataVT, aes(log(V), log(Thickness), color=Suborder)) +
geom_point(size=4) +
xlab("ln egg volume (mm3)") +
ylab("ln eggshell thickness (µm)") +
geom_abline(intercept=NewEBVT$coefficients[1], slope=NewEBVT$coefficients[2], colour="lightblue") +
theme(panel.background = element_rect(fill="black")) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_colour_brewer("Suborder", palette="Accent")
ppadata<-dataVT
ppadata[,c(4,6:8)]<-log(ppadata[,c(4,6:8)])
M<-define_model_set(null=c(),
one=c(Thickness~SVL),
two=c(V~BM),
three=c(Thickness~V),
four=c(Thickness~V+SVL),
five=c(Thickness~V, V~BM),
six=c(Thickness~SVL, V~BM),
seven=c(Thickness~SVL+BM),
eight=c(V~BM, Thickness~BM),
nine=c(Thickness~V+BM),
ten=c(Thickness~SVL+BM+V),
eleven=c(Thickness~V+BM, V~BM),
twelve=c(Thickness~SVL+BM, V~BM),
thirteen=c(Thickness~BM),
.common=c(BM~SVL, V~SVL))
plot_model_set(M)
## Warning: The curvature argument has been deprecated in favour of strength
result<-phylo_path(M, data=ppadata, tree=pruned.treeVT, model='lambda'); result
## A phylogenetic path analysis, on the variables:
## Continuous: Thickness SVL V BM
## Binary:
##
## Evaluated for these models: null one two three four five six seven eight nine ten eleven twelve thirteen
##
## Containing 31 phylogenetic regressions, of which 10 unique
s<-summary(result); s
plot(s)
## Warning: Position guide is perpendicular to the intended axis. Did you mean to
## specify a different guide `position`?
averagemodel<-average(result)
## Registered S3 methods overwritten by 'MuMIn':
## method from
## nobs.pgls caper
## nobs.phylolm phylolm
## logLik.phylolm phylolm
plot(averagemodel)
## Warning: The curvature argument has been deprecated in favour of strength
fossiltreePEM<-read.nexus("Lepidosaurtree.trees.nex")
fossildataPEM<-read.table("Dataset1-lepidosaurs.txt",header=TRUE,stringsAsFactor=FALSE,)
fossildataV<-subset(fossildataPEM, !is.na(V))
fossildataV[,c(3:5,7:9)]<-log(fossildataV[,c(3:5,7:9)])
treePEM<-drop.tip(fossiltreePEM, setdiff(fossiltreePEM$tip.label, fossildataV$Species))
treedrop <- drop.tip(treePEM,fossildataV[is.na(fossildataV[,"SVL"]),"Species"])
grloc <- getGraphLocations(treePEM,fossildataV[is.na(fossildataV[,"SVL"]),"Species"])
sporder <- match(attr(grloc$x,"vlabel")[grloc$x$vertex$species],fossildataV[,"Species"])
PEMfs <- list()
PEMfs[["V"]] <- PEM.fitSimple(y=fossildataV[sporder,"SVL"],
x=fossildataV[sporder,"V"],w=grloc$x,d="distance",sp="species",
lower=0,upper=1)
PEMfs[["none"]] <- PEM.fitSimple(y=fossildataV[sporder,"SVL"],
x=NULL,w=grloc$x,d="distance",sp="species",lower=0,upper=1)
for(m in c("V","none")) print(PEMfs[[m]]$optim$par)
## [1] 0.4981507
## [1] 0.1197493
rm(m)
PEMAIC <- list()
PEMAIC[["V"]] <- lmforwardsequentialAICc(y=fossildataV[sporder,"SVL"],
x=fossildataV[sporder,"V",drop=FALSE],object=PEMfs[["V"]])
PEMAIC[["none"]] <- lmforwardsequentialAICc(y=fossildataV[sporder,"SVL"],object=PEMfs[["none"]])
for(m in c("V","none"))
cat(m,summary(PEMAIC[[m]])$adj,PEMAIC[[m]]$AICc,"\n")
## V 0.9945104 -248.9384
## none 0.9913577 -124.1351
rm(m)
summary(PEMAIC[["V"]])
##
## Call:
## lm(formula = as.formula(paste(p1, p2, sep = "")), data = df1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.178636 -0.035423 0.004214 0.039172 0.153809
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.071635 0.066712 31.053 < 2e-16 ***
## V 0.393493 0.009214 42.705 < 2e-16 ***
## V_1 -9.100852 0.108193 -84.116 < 2e-16 ***
## V_2 -4.704946 0.126827 -37.097 < 2e-16 ***
## V_30 1.372985 0.083670 16.410 < 2e-16 ***
## V_6 0.986086 0.110409 8.931 3.32e-15 ***
## V_10 0.858355 0.118694 7.232 3.56e-11 ***
## V_8 -0.984514 0.083734 -11.758 < 2e-16 ***
## V_63 -0.760769 0.084320 -9.022 1.98e-15 ***
## V_7 0.862697 0.083710 10.306 < 2e-16 ***
## V_31 0.874262 0.084017 10.406 < 2e-16 ***
## V_117 0.819127 0.084435 9.701 < 2e-16 ***
## V_27 -1.240447 0.106374 -11.661 < 2e-16 ***
## V_77 -0.733867 0.083703 -8.768 8.30e-15 ***
## V_229 -0.819382 0.086023 -9.525 < 2e-16 ***
## V_156 0.507389 0.085388 5.942 2.39e-08 ***
## V_124 -0.408268 0.089180 -4.578 1.08e-05 ***
## V_5 0.388503 0.087559 4.437 1.91e-05 ***
## V_180 -0.572372 0.083671 -6.841 2.72e-10 ***
## V_14 -0.450315 0.085040 -5.295 4.87e-07 ***
## V_225 -0.546724 0.083745 -6.528 1.33e-09 ***
## V_48 -0.617529 0.083862 -7.364 1.77e-11 ***
## V_11 0.632745 0.084051 7.528 7.39e-12 ***
## V_108 0.675610 0.084757 7.971 6.76e-13 ***
## V_76 0.648948 0.084352 7.693 3.05e-12 ***
## V_44 -0.384803 0.085978 -4.476 1.64e-05 ***
## V_17 0.485665 0.083976 5.783 5.10e-08 ***
## V_80 0.535035 0.083675 6.394 2.61e-09 ***
## V_205 0.551436 0.083744 6.585 1.00e-09 ***
## V_23 -0.401879 0.084055 -4.781 4.61e-06 ***
## V_59 -0.320047 0.085352 -3.750 0.000265 ***
## V_66 0.544710 0.084217 6.468 1.81e-09 ***
## V_13 -0.499332 0.083867 -5.954 2.26e-08 ***
## V_29 -0.594629 0.085593 -6.947 1.57e-10 ***
## V_170 0.403913 0.083791 4.820 3.90e-06 ***
## V_109 -0.541708 0.084637 -6.400 2.54e-09 ***
## V_98 -0.558359 0.085053 -6.565 1.11e-09 ***
## V_179 -0.390290 0.083857 -4.654 7.86e-06 ***
## V_18 -0.529632 0.085127 -6.222 6.15e-09 ***
## V_206 -0.439153 0.083739 -5.244 6.13e-07 ***
## V_54 -0.453173 0.083883 -5.402 3.00e-07 ***
## V_125 -0.458915 0.084018 -5.462 2.28e-07 ***
## V_142 0.402447 0.083669 4.810 4.08e-06 ***
## V_235 0.429425 0.083804 5.124 1.05e-06 ***
## V_155 0.370097 0.083716 4.421 2.04e-05 ***
## V_34 0.351834 0.083814 4.198 4.94e-05 ***
## V_194 0.273009 0.085159 3.206 0.001692 **
## V_183 -0.376601 0.083670 -4.501 1.48e-05 ***
## V_169 -0.344993 0.083766 -4.119 6.71e-05 ***
## V_112 0.278725 0.084640 3.293 0.001274 **
## V_39 -0.293362 0.084304 -3.480 0.000682 ***
## V_55 -0.391209 0.083709 -4.673 7.26e-06 ***
## V_69 0.380311 0.083680 4.545 1.24e-05 ***
## V_129 -0.395590 0.083754 -4.723 5.89e-06 ***
## V_173 -0.381593 0.083692 -4.559 1.16e-05 ***
## V_152 -0.131500 0.088366 -1.488 0.139119
## V_56 0.147425 0.087294 1.689 0.093629 .
## V_62 0.424582 0.084164 5.045 1.48e-06 ***
## V_104 0.367887 0.083752 4.393 2.29e-05 ***
## V_72 -0.370855 0.083789 -4.426 2.00e-05 ***
## V_88 -0.321197 0.083667 -3.839 0.000191 ***
## V_178 0.381534 0.084033 4.540 1.26e-05 ***
## V_133 -0.435337 0.084932 -5.126 1.04e-06 ***
## V_58 0.361402 0.083871 4.309 3.19e-05 ***
## V_68 0.363490 0.084099 4.322 3.03e-05 ***
## V_172 -0.275813 0.083706 -3.295 0.001266 **
## V_176 0.264259 0.083676 3.158 0.001972 **
## V_151 0.262899 0.083679 3.142 0.002077 **
## V_101 -0.279064 0.083671 -3.335 0.001109 **
## V_130 0.354518 0.084439 4.199 4.93e-05 ***
## V_20 0.450383 0.087594 5.142 9.68e-07 ***
## V_16 0.364896 0.084960 4.295 3.38e-05 ***
## V_4 0.448939 0.088214 5.089 1.22e-06 ***
## V_134 -0.322932 0.084261 -3.833 0.000196 ***
## V_45 -0.265568 0.083671 -3.174 0.001874 **
## V_120 0.273162 0.083703 3.263 0.001404 **
## V_227 -0.273278 0.083759 -3.263 0.001408 **
## V_33 0.347907 0.085649 4.062 8.33e-05 ***
## V_148 0.267991 0.083752 3.200 0.001725 **
## V_118 -0.283059 0.083948 -3.372 0.000982 ***
## V_138 0.267968 0.083770 3.199 0.001731 **
## V_121 -0.253752 0.083678 -3.032 0.002925 **
## V_42 0.299744 0.084538 3.546 0.000544 ***
## V_97 -0.272391 0.083932 -3.245 0.001489 **
## V_116 -0.260448 0.083813 -3.107 0.002314 **
## V_73 -0.209929 0.083934 -2.501 0.013613 *
## V_85 0.228299 0.083680 2.728 0.007242 **
## V_19 0.264901 0.084013 3.153 0.002003 **
## V_28 -0.307979 0.086077 -3.578 0.000486 ***
## V_26 0.224630 0.083682 2.684 0.008207 **
## V_21 -0.292578 0.085859 -3.408 0.000871 ***
## V_61 -0.242060 0.083830 -2.888 0.004544 **
## V_110 -0.242602 0.083878 -2.892 0.004479 **
## V_177 0.245034 0.084113 2.913 0.004208 **
## V_193 -0.215603 0.083688 -2.576 0.011095 *
## V_99 -0.244967 0.085178 -2.876 0.004704 **
## V_95 -0.219343 0.084007 -2.611 0.010080 *
## V_187 0.201047 0.083676 2.403 0.017677 *
## V_115 -0.194921 0.083706 -2.329 0.021409 *
## V_40 -0.204034 0.083748 -2.436 0.016183 *
## V_35 0.231018 0.086910 2.658 0.008837 **
## V_38 0.223027 0.085936 2.595 0.010530 *
## V_150 -0.203677 0.084038 -2.424 0.016732 *
## V_222 -0.198031 0.083780 -2.364 0.019564 *
## V_168 -0.188724 0.083980 -2.247 0.026296 *
## V_141 0.183297 0.083777 2.188 0.030449 *
## V_12 0.180644 0.083706 2.158 0.032745 *
## V_91 -0.174506 0.084069 -2.076 0.039873 *
## V_201 -0.168706 0.083691 -2.016 0.045863 *
## V_67 -0.162550 0.083667 -1.943 0.054182 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08367 on 131 degrees of freedom
## Multiple R-squared: 0.997, Adjusted R-squared: 0.9945
## F-statistic: 399.9 on 109 and 131 DF, p-value: < 2.2e-16
The model with V as a co-predictor has the best fit.
tf<-list(V=identity, none=identity)
m <- "V" ; atr <- tf[[m]](fossildataV[is.na(fossildataV[,"SVL"]),m,drop=FALSE])
resultsPEM<-predict(object=PEMfs[[m]],targets=grloc,lmobject=PEMAIC[[m]],newdata=atr,interval="confidence")
# Predicted SVL, with upper and lower limits of the confidence interval for the prediction:
exp(resultsPEM$values); exp(resultsPEM$upper); exp(resultsPEM$lower)
## [,1]
## Antarcticoolithus_bradyi 6665.867
## [,1]
## Antarcticoolithus_bradyi 7888.606
## [,1]
## Antarcticoolithus_bradyi 5632.654
The code to predict missing values for body mass is identical to the one for SVL; simply replace ‘SVL’ with ‘BM’.
data[1,7]<-6665.867
dataV2<-subset(data, !is.na(V))
ggplot(dataV2, aes(log(SVL), log(V), color=Suborder)) +
geom_point(size=4) +
xlab("ln snout-vent length (mm)") +
ylab("ln egg volume (mm3)") +
geom_abline(intercept=Lambda$coefficients[1], slope=Lambda$coefficients[2], colour="lightblue") +
theme(panel.background = element_rect(fill="black")) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_colour_brewer("Suborder", palette="Accent")
data[1,8]<-456096.7
dataV2<-subset(data, !is.na(V))
ggplot(dataV2, aes(log(BM), log(V), color=Suborder)) +
geom_point(size=4) +
xlab("ln body mass (g)") +
ylab("ln egg volume (mm3)") +
geom_abline(intercept=LambdaBM$coefficients[1], slope=LambdaBM$coefficients[2], colour="lightblue") +
theme(panel.background = element_rect(fill="black")) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_colour_brewer("Suborder", palette="Set1")
RT<-data$Thickness/data$L; data<-cbind(data, RT)
dataA<-subset(data, !is.na(RT)&!is.na(OV))
dataA$RT<-log(dataA$RT)
anovatree<-drop.tip(tree, setdiff(tree$tip.label, rownames(dataA)))
datANOVA<-as.data.frame(cbind(dataA$RT, dataA$OV), header=T)
row.names(datANOVA)<-row.names(dataA); colnames(datANOVA)<-c("RT", "OV")
datANOVA<-ReorderData(anovatree, datANOVA, taxa.names="row names")
phylANOVA(anovatree, datANOVA$OV, datANOVA$RT, p.adj="BH")
## Warning: no labels for x. Assuming order of tree$tip.label.
##
## Warning: no labels for y. Assuming order of tree$tip.label.
## ANOVA table: Phylogenetic ANOVA
##
## Response: y
## Sum Sq Mean Sq F value Pr(>F)
## x 2.884385 2.884385 6.074072 0.011
## Residual 31.816186 0.474868
##
## P-value based on simulation.
## ---------
##
## Pairwise posthoc test using method = "BH"
##
## Pairwise t-values:
## 1 2
## 1 0.000000 2.464563
## 2 -2.464563 0.000000
##
## Pairwise corrected P-values:
## 1 2
## 1 1.000 0.011
## 2 0.011 1.000
## ---------
ggplot(dataA, aes(x=OV, y=RT, fill=OV)) +
geom_boxplot() +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_fill_hue()
log(data$RT)[1] # Value for Antarcticoolithus, intermediate between the groups
## [1] 0.9120565
p2<-ggplot(dataA, aes(x=OV, y=RT, fill=OV)) + geom_boxplot()
ggplot_build(p2)$data
## [[1]]
## fill ymin lower middle upper ymax outliers
## 1 #F8766D -0.1251631 0.9808293 1.35720692 1.7429693 2.784423 3.442019
## 2 #00BFC4 -0.3462762 -0.3308263 -0.02868824 0.8270896 2.503459
## notchupper notchlower x flipped_aes PANEL group ymin_final ymax_final xmin
## 1 1.5065672 1.2078466 1 FALSE 1 1 -0.1251631 3.442019 0.625
## 2 0.8860653 -0.9434418 2 FALSE 1 2 -0.3462762 2.503459 1.625
## xmax xid newx new_width weight colour size alpha shape linetype
## 1 1.375 1 1 0.75 1 grey20 0.5 NA 19 solid
## 2 2.375 2 2 0.75 1 grey20 0.5 NA 19 solid
datANOVA2<-as.data.frame(cbind(dataA$RT, dataA$Strategy), header=T)
row.names(datANOVA2)<-row.names(dataA); colnames(datANOVA2)<-c("RT", "Strategy")
datANOVA2<-ReorderData(anovatree, datANOVA2, taxa.names="row names")
phylANOVA(anovatree, datANOVA2$Strategy, datANOVA2$RT, p.adj="BH")
## Warning: no labels for x. Assuming order of tree$tip.label.
##
## Warning: no labels for y. Assuming order of tree$tip.label.
## ANOVA table: Phylogenetic ANOVA
##
## Response: y
## Sum Sq Mean Sq F value Pr(>F)
## x 7.86154 3.930770 9.666176 0.005
## Residual 26.83903 0.406652
##
## P-value based on simulation.
## ---------
##
## Pairwise posthoc test using method = "BH"
##
## Pairwise t-values:
## 1 2 3
## 1 0.000000 -1.534310 3.960211
## 2 1.534310 0.000000 4.382854
## 3 -3.960211 -4.382854 0.000000
##
## Pairwise corrected P-values:
## 1 2 3
## 1 1.000 0.311 0.003
## 2 0.311 1.000 0.003
## 3 0.003 0.003 1.000
## ---------
ggplot(dataA, aes(x=Strategy, y=RT, fill=Strategy)) +
geom_boxplot() +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_fill_hue()
p2<-ggplot(dataA, aes(x=Strategy, y=RT, fill=Strategy)) + geom_boxplot()
ggplot_build(p2)$data
## [[1]]
## fill ymin lower middle upper ymax outliers
## 1 #F8766D -0.1251631 0.9794977 1.3139737 1.74296931 2.7844232
## 2 #00BA38 0.3137979 1.3231783 1.6757606 2.06667501 2.5034588 3.442019
## 3 #619CFF -0.3462762 -0.3359763 -0.3256763 -0.02868824 0.2682998
## notchupper notchlower x flipped_aes PANEL group ymin_final ymax_final xmin
## 1 1.47662916 1.1513182 1 FALSE 1 1 -0.1251631 2.7844232 0.625
## 2 2.02995340 1.3215677 2 FALSE 1 2 0.3137979 3.4420194 1.625
## 3 -0.04536403 -0.6059886 3 FALSE 1 3 -0.3462762 0.2682998 2.625
## xmax xid newx new_width weight colour size alpha shape linetype
## 1 1.375 1 1 0.75 1 grey20 0.5 NA 19 solid
## 2 2.375 2 2 0.75 1 grey20 0.5 NA 19 solid
## 3 3.375 3 3 0.75 1 grey20 0.5 NA 19 solid
This part requires the use of dataset 2 (“Dataset2-amniotes.txt”; see also Supplementary Table 3 in Legendre et al., 2020), and amniote phylogenetic tree (“Amniotetree.nex”).
dataS<-read.table("Dataset2-amniotes.txt", header=T)
treeS<-read.nexus("Amniotetree.nex")
treeS<-drop.tip(treeS, setdiff(treeS$tip.label, dataS$Taxon))
dataS<-read.table("Dataset2-amniotes.txt", header=T, row.names="Taxon")
Ws<-diag(vcv.phylo(treeS))
dataS<-read.table("Dataset2-amniotes.txt", header=T)
datahard<-dataS[-c(50,86:88,93:95,115,118,120:145,147,149),]
treehard<-drop.tip(treeS, setdiff(treeS$tip.label, datahard$Taxon))
dataS<-read.table("Dataset2-amniotes.txt", header=T, row.names="Taxon")
datahard<-dataS[-c(50,86:88,93:95,115,118,120:145,147,149),]
Wh<-diag(vcv.phylo(treehard))
alpha <- seq(0, 1, 0.1)
fit <- list()
form <- log(Eggshell_thickness)~log(Egg_mass)
for (i in seq_along(alpha)) {
cor <- corMartins(alpha[i], phy = treehard, fixed = T)
fit[[i]] <- gls(form, correlation = cor, data = datahard, na.action=na.exclude, weights=varFixed(~Wh), method = "ML")
}
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
sapply(fit, logLik)
## [1] Inf -37.10514 -37.00214 -37.38638 -37.59015 -37.68151 -37.71941
## [8] -37.73459 -37.74057 -37.74291 -37.74382
g <- seq(0.1, 1, 0.1)
fit <- list()
form <- log(Eggshell_thickness)~log(Egg_mass)
for (i in seq_along(g)) {
cor <- corBlomberg(g[i], phy = treehard, fixed = T)
fit[[i]] <- gls(form, correlation = cor, data = datahard, na.action=na.exclude, weights=varFixed(~Wh), method = "ML")
}
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
sapply(fit, logLik)
## [1] -51.97396 -62.61195 -68.20028 -71.62168 -73.95198 -75.65088 -76.94873
## [8] -77.97487 -78.80802 -79.49898
BM<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datahard, correlation=corBrownian(phy=treehard), weights=varFixed(~Wh), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OU<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datahard, correlation=corMartins(0.1, phy=treehard, fixed=T), weights=varFixed(~Wh), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
Lambda<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datahard, correlation=corPagel(1, phy=treehard), weights=varFixed(~Wh), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
EB<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datahard, correlation=corBlomberg(0.1, phy=treehard, fixed=T), weights=varFixed(~Wh), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OLS<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datahard, method="ML")
Cand.models = list()
Cand.models[[1]] = BM
Cand.models[[2]] = OU
Cand.models[[3]] = Lambda
Cand.models[[4]] = EB
Cand.models[[5]] = OLS
Modnames = paste(c("BM", "OU", "Lambda", "EB", "OLS"), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = T)
summary(Lambda)
## Generalized least squares fit by maximum likelihood
## Model: log(Eggshell_thickness) ~ log(Egg_mass)
## Data: datahard
## AIC BIC logLik
## 67.05542 77.92942 -29.52771
##
## Correlation Structure: corPagel
## Formula: ~1
## Parameter estimate(s):
## lambda
## 0.3659904
## Variance function:
## Structure: fixed weights
## Formula: ~Wh
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 4.386080 0.12772437 34.34020 0
## log(Egg_mass) 0.426746 0.01576025 27.07737 0
##
## Correlation:
## (Intr)
## log(Egg_mass) -0.443
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.78157887 -0.99474516 -0.57993810 0.03130181 2.28401650
##
## Residual standard error: 0.02035964
## Degrees of freedom: 112 total; 110 residual
caperdataS<-read.table("Dataset2-amniotes.txt", header=T)
datahard<-dataS[-c(50,86:88,93:95,115,118,120:145,147,149),]
datacomp<-comparative.data(phy=treehard, data=datahard, names.col="Taxon")
pgls<-pgls(log(Eggshell_thickness)~log(Egg_mass), data=datacomp, lambda="ML")
summary(pgls)
##
## Call:
## pgls(formula = log(Eggshell_thickness) ~ log(Egg_mass), data = datacomp,
## lambda = "ML")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.06443 -0.01539 -0.00283 0.01159 0.05560
##
## Branch length transformations:
##
## kappa [Fix] : 1.000
## lambda [ ML] : 0.328
## lower bound : 0.000, p = 0.0012278
## upper bound : 1.000, p = < 2.22e-16
## 95.0% CI : (0.098, 0.597)
## delta [Fix] : 1.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.334504 0.128404 33.757 < 2.2e-16 ***
## log(Egg_mass) 0.436835 0.017409 25.092 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.02103 on 110 degrees of freedom
## Multiple R-squared: 0.8513, Adjusted R-squared: 0.8499
## F-statistic: 629.6 on 1 and 110 DF, p-value: < 2.2e-16
datasoft<-dataS[c(50,86:88,93,94,115,118,120:145,147,149),]
treesoft<-drop.tip(treeS, setdiff(treeS$tip.label, datasoft$Taxon))
dataS<-read.table("Dataset2-amniotes.txt", header=T, row.names="Taxon")
datasoft<-dataS[c(50,86:88,93,94,115,118,120:145,147,149),]
Ws<-diag(vcv.phylo(treesoft))
alpha <- seq(0, 1, 0.1)
fit <- list()
form <- log(Eggshell_thickness)~log(Egg_mass)
for (i in seq_along(alpha)) {
cor <- corMartins(alpha[i], phy = treesoft, fixed = T)
fit[[i]] <- gls(form, correlation = cor, data = datasoft, na.action=na.exclude, weights=varFixed(~Ws), method = "ML")
}
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
sapply(fit, logLik)
## [1] Inf -31.45442 -31.58754 -31.70751 -31.76148 -31.78240 -31.79013
## [8] -31.79293 -31.79394 -31.79430 -31.79443
g <- seq(0.1, 1, 0.1)
fit <- list()
form <- log(Eggshell_thickness)~log(Egg_mass)
for (i in seq_along(g)) {
cor <- corBlomberg(g[i], phy = treesoft, fixed = T)
fit[[i]] <- gls(form, correlation = cor, data = datasoft, na.action=na.exclude, weights=varFixed(~Ws), method = "ML")
}
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
sapply(fit, logLik)
## [1] -31.02830 -33.27681 -35.31325 -36.86832 -38.06022 -38.99694 -39.75253
## [8] -40.37614 -40.90094 -41.34988
BM<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datasoft, correlation=corBrownian(phy=treesoft), weights=varFixed(~Ws), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OU<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datasoft, correlation=corMartins(0.2, phy=treesoft), weights=varFixed(~Ws), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
## Error in recalc.corStruct(object[[i]], conLin): NA/NaN/Inf in foreign function call (arg 4)
Lambda<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datasoft, correlation=corPagel(1, phy=treesoft), weights=varFixed(~Ws), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
EB<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datasoft, correlation=corBlomberg(0.1, phy=treesoft, fixed=T), weights=varFixed(~Ws), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
OLS<-gls(log(Eggshell_thickness)~log(Egg_mass), data=datasoft, method="ML")
Cand.models = list()
Cand.models[[1]] = BM
Cand.models[[2]] = OU
Cand.models[[3]] = Lambda
Cand.models[[4]] = EB
Cand.models[[5]] = OLS
Modnames = paste(c("BM", "OU", "Lambda", "EB", "OLS"), sep = " ")
aictab(cand.set = Cand.models, modnames = Modnames, sort = T)
summary(EB)
## Generalized least squares fit by maximum likelihood
## Model: log(Eggshell_thickness) ~ log(Egg_mass)
## Data: datasoft
## AIC BIC logLik
## 68.05659 72.80715 -31.0283
##
## Correlation Structure: corBlomberg
## Formula: ~1
## Parameter estimate(s):
## g
## 0.1
## Variance function:
## Structure: fixed weights
## Formula: ~Ws
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 2.0167850 0.16983372 11.875056 0
## log(Egg_mass) 0.4009087 0.04963276 8.077501 0
##
## Correlation:
## (Intr)
## log(Egg_mass) -0.643
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.0217026 -0.6275853 -0.1318103 0.6864130 1.6585487
##
## Residual standard error: 0.03771851
## Degrees of freedom: 36 total; 34 residual
# R-squared
EB2<-gls(log(Eggshell_thickness)~1, data=datasoft, correlation=corBlomberg(0.1, phy=treesoft, fixed=T), weights=varFixed(~Ws), method="ML")
## Warning in Initialize.corPhyl(X[[i]], ...): No covariate specified, species
## will be taken as ordered in the data frame. To avoid this message, specify a
## covariate containing the species names with the 'form' argument.
1 - (EB$sigma/EB2$sigma)^2
## [1] 0.657417
# P-value
anova(EB2, EB)
dataS<-read.table("Dataset2-amniotes.txt", header=T)
ggplot(dataS, aes(log(Egg_mass), log(Eggshell_thickness), colour=Clade)) +
geom_point(size=4) +
xlab("ln egg mass (g)") +
ylab("ln calcareous layer thickness (µm)") +
geom_abline(intercept=4.321066, slope=0.439690, colour="turquoise", size=1.3) +
geom_abline(intercept=2.0474668, slope=0.4036309, colour="red", size=1.3) +
theme(panel.background = element_rect(fill="black")) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_colour_brewer("Group", palette="Set1")
fastAncdataS<-read.table("Dataset2-amniotes.txt", header=T); dataS<-dataS[-95,]
treeplot<-drop.tip(treeS, setdiff(treeS$tip.label, dataS$Taxon))
dataS<-read.table("Dataset2-amniotes.txt", header=T, row.names="Taxon"); dataS<-dataS[-95,]
dataplot=as.matrix(dataS[,3]); names(dataplot)<-rownames(dataS); colnames(dataplot)<-"CL_thickness"
dataplot<-log10(dataplot)
fit<-fastAnc(treeplot,dataplot,vars=T,CI=T) # Ancestral states for each node
obj<-contMap(treeplot,dataplot,plot=F)
plot(setMap(obj, colors=rev(brewer.pal(10, "Spectral"))), fsize=0.4,lwd=4)
treeplot<-force.ultrametric(treeplot)
fit<-fastAnc(treeplot,dataplot,vars=T,CI=T) # Ancestral states for each node
obj<-contMap(treeplot,dataplot,plot=F)
plot(setMap(obj, colors=rev(brewer.pal(10, "Spectral"))), fsize=0.4,lwd=4)
make.simmapdataS<-read.table("Dataset2-amniotes.txt", header=T)
dataS<-dataS[-95,]
treeS<-read.nexus("Amniotetree.nex")
treeS<-drop.tip(treeS, setdiff(treeS$tip.label, dataS$Taxon))
prismatic<-dataS[,5]; names(prismatic)<-dataS$Taxon
cols<-setNames(c("royalblue","red3"),levels(prismatic))
mtree<-make.simmap(treeS, prismatic, model="ER")
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## Absent Present
## Absent -0.0003949281 0.0003949281
## Present 0.0003949281 -0.0003949281
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Absent Present
## 0.5 0.5
## Done.
plot(mtree,cols,fsize=0.5,ftype="i")
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
y=0.8*par()$usr[3],fsize=0.8)
# with 1000 simulations of stochastic character maps from the data:
mtrees<-make.simmap(treeS, prismatic, model="ER", nsim=1000)
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## Absent Present
## Absent -0.0003949281 0.0003949281
## Present 0.0003949281 -0.0003949281
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## Absent Present
## 0.5 0.5
## Done.
par(mfrow=c(10,10))
null<-sapply(mtrees,plotSimmap,colors=cols,lwd=1,ftype="off")
pd<-summary(mtrees); pd
## 1000 trees with a mapped discrete character with states:
## Absent, Present
##
## trees have 3.735 changes between states on average
##
## changes are of the following types:
## Absent,Present Present,Absent
## x->y 0.981 2.754
##
## mean total time spent in each state is:
## Absent Present total
## raw 2705.1756568 6756.0840379 9461.26
## prop 0.2859213 0.7140787 1.00
par(mfrow=c(1,1))
plot(pd,fsize=0.6,ftype="i",colors=cols,ylim=c(-2,Ntip(treeS)))
add.simmap.legend(colors=cols[2:1],prompt=FALSE,x=0,y=-2,vertical=FALSE)
# density map of the above
objprism<-densityMap(mtrees,states=levels(prismatic)[2:1],plot=FALSE)
## sorry - this might take a while; please be patient
n<-length(objprism$cols)
objprism$cols[1:n]<-colorRampPalette(c("royalblue","red3"))(n)
plot(objprism,size=c(0.6,1),fsize=0.6,lwd=4)
dataS<-read.table("Dataset2-amniotes.txt", header=T)
RTM<-dataS$Eggshell_thickness/dataS$Egg_mass
dataS<-cbind(dataS, RTM)
dataSboxplot<-dataS %>% filter(!Clade=="Test")
ggplot(dataSboxplot, aes(x=Clade, y=log(RTM), fill=Clade)) +
geom_boxplot() +
theme(panel.background = element_rect(fill="white")) +
theme(axis.text=element_text(size=10, angle = 90),
axis.title=element_text(size=15),
legend.text=element_text(size=15),
legend.title=element_text(size=17)) +
scale_fill_brewer(palette="Set1")
p<-ggplot(dataSboxplot, aes(x=Clade, y=log(RTM), fill=Clade)) +
geom_boxplot()
ggplot_build(p)$data
## [[1]]
## fill ymin lower middle upper ymax
## 1 #F8766D -4.1737717 -4.1737717 -4.1737717 -4.1737717 -4.1737717
## 2 #D39200 0.4296871 1.8265428 2.4120341 2.9394662 4.1207778
## 3 #93AA00 0.2636401 1.6216397 2.5162801 2.5794171 2.9038769
## 4 #00BA38 1.4714471 1.4714471 1.4714471 1.4714471 1.4714471
## 5 #00C19F 1.1282254 1.5449861 1.7324607 1.9034789 2.2137481
## 6 #00B9E3 4.4770380 4.4770380 4.4770380 4.4770380 4.4770380
## 7 #619CFF -2.1102132 -0.4793495 -0.1364957 0.9978001 1.9902104
## 8 #DB72FB 0.2876821 0.3729137 0.4581454 0.5433770 0.6286087
## 9 #FF61C3 -1.4307461 0.2392498 1.7701329 3.4510661 5.4647026
## outliers notchupper notchlower x flipped_aes y PANEL
## 1 -4.1737717 -4.1737717 1 FALSE -4.173772 1
## 2 -0.9350988, -0.6852781 2.6607121 2.1633561 2 FALSE NA 1
## 3 3.0207095 2.0118506 3 FALSE NA 1
## 4 1.4714471 1.4714471 4 FALSE 1.471447 1
## 5 1.9327199 1.5322015 5 FALSE NA 1
## 6 4.4770380 4.4770380 6 FALSE 4.477038 1
## 7 0.2372268 -0.5102183 7 FALSE NA 1
## 8 0.6485919 0.2676989 8 FALSE NA 1
## 9 2.5933530 0.9469128 9 FALSE NA 1
## group ymin_final ymax_final xmin xmax xid newx new_width weight colour size
## 1 1 -4.1737717 -4.1737717 0.625 1.375 1 1 0.75 1 grey20 0.5
## 2 2 -0.9350988 4.1207778 1.625 2.375 2 2 0.75 1 grey20 0.5
## 3 3 0.2636401 2.9038769 2.625 3.375 3 3 0.75 1 grey20 0.5
## 4 4 1.4714471 1.4714471 3.625 4.375 4 4 0.75 1 grey20 0.5
## 5 5 1.1282254 2.2137481 4.625 5.375 5 5 0.75 1 grey20 0.5
## 6 6 4.4770380 4.4770380 5.625 6.375 6 6 0.75 1 grey20 0.5
## 7 7 -2.1102132 1.9902104 6.625 7.375 7 7 0.75 1 grey20 0.5
## 8 8 0.2876821 0.6286087 7.625 8.375 8 8 0.75 1 grey20 0.5
## 9 9 -1.4307461 5.4647026 8.625 9.375 9 9 0.75 1 grey20 0.5
## alpha shape linetype
## 1 NA 19 solid
## 2 NA 19 solid
## 3 NA 19 solid
## 4 NA 19 solid
## 5 NA 19 solid
## 6 NA 19 solid
## 7 NA 19 solid
## 8 NA 19 solid
## 9 NA 19 solid